Bayesian GLM Part1

Preparations

Load the necessary libraries

library(rstanarm)   #for fitting models in STAN
library(brms)       #for fitting models in STAN
library(standist)   #for exploring distributions
library(coda)       #for diagnostics
library(bayesplot)  #for diagnostics
library(ggmcmc)     #for MCMC diagnostics
library(DHARMa)     #for residual diagnostics
library(rstan)      #for interfacing with STAN
library(emmeans)    #for marginal means etc
library(broom)      #for tidying outputs
library(tidybayes)  #for more tidying outputs
library(ggeffects)  #for partial plots
library(tidyverse)  #for data wrangling etc
library(broom.mixed)#for summarising models
library(ggeffects)  #for partial effects plots
theme_set(theme_grey()) #put the default ggplot theme back

Scenario

Here is an example from Fowler, Cohen, and Jarvis (1998). An agriculturalist was interested in the effects of fertilizer load on the yield of grass. Grass seed was sown uniformly over an area and different quantities of commercial fertilizer were applied to each of ten 1 m² randomly located plots. Two months later the grass from each plot was harvested, dried and weighed. The data are in the file fertilizer.csv in the data folder.

FERTILIZER	YIELD
25	84
50	80
75	90
100	154
125	148
...	...

FERTILIZER:	Mass of fertilizer (g.m-2) - Predictor variable
YIELD:	Yield of grass (g.m-2) - Response variable

The aim of the analysis is to investigate the relationship between fertilizer concentration and grass yield.

Read in the data

We will start off by reading in the Fertilizer data. There are many functions in R that can read in a CSV file. We will use a the read_csv() function as it is part of the tidyverse ecosystem.

fert = read_csv('../public/data/fertilizer.csv', trim_ws=TRUE)
glimpse(fert)

## Rows: 10
## Columns: 2
## $ FERTILIZER <dbl> 25, 50, 75, 100, 125, 150, 175, 200, 225, 250
## $ YIELD      <dbl> 84, 80, 90, 154, 148, 169, 206, 244, 212, 248

Exploratory data analysis

The individual responses (\(y_i\), observed yields) are each expected to have been independently drawn from normal (Gaussian) distributions (\(\mathcal{N}\)). These distributions represent all the possible yields we could have obtained at the specific (\(i^th\)) level of Fertilizer. Hence the \(i^th\) yield observation is expected to have been drawn from a normal distribution with a mean of \(\mu_i\).

Although each distribution is expected to come from populations that differ in their means, we assume that all of these distributions have the same variance (\(\sigma^2\)).

We need to supply priors for each of the parameters to be estimated (\(\beta_0\), \(\beta_1\) and \(\sigma\)). Whilst we want these priors to be sufficiently vague as to not influence the outcomes of the analysis (and thus be equivalent to the frequentist analysis), we do not want the priors to be so vague (wide) that they permit the MCMC sampler to drift off into parameter space that is both ilogical as well as numerically awkward.

As a starting point, lets assign the following priors:

• \(\beta_0\): Normal prior centered at 0 with a variance of 100
• \(\beta_1\): Normal prior centered at 0 with a variance of 10
• \(\sigma\): either a half cauchy centered at 0 with a variance of 5 or an exponential with parameter of 1

Note, when fitting models through either rstanarm or brms, the priors assume that the predictor(s) have been centered and are to be applied on the link scale. In this case the link scale is an identity.

Model formula: \[ \begin{align} y_i &\sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i &= \beta_0 + \beta_1 x_i\\ \beta_0 &\sim{} \mathcal{N}(0,100)\\ \beta_1 &\sim{} \mathcal{N}(0,10)\\ \sigma &\sim{} \mathcal{cauchy}(0,5)\\ OR\\ \sigma &\sim{} \mathcal{Exp}(1)\\ \end{align} \]

Fit the model

Note, for routines that take more than a couple of seconds to perform (such as most Bayesian models), it is a good idea to cache the Rmarkdown chunks in which the routine is performed. That way, the routine will only be run the first time and any objects generated will be stored for future use. Thereafter, provided the code has not changed, the routine will not be re-run. Rather, knitr will just retrieve the cached objects and continue on.

frequentist

For the purpose of comparing the frequentist and Bayesian outcomes, it might be useful to start by fitting the frequestist simple linear model. We will also fit this model with the predictor (fertilizer concentration) centered.

summary(lm(YIELD ~ FERTILIZER, data=fert))

## 
## Call:
## lm(formula = YIELD ~ FERTILIZER, data = fert)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -22.79 -11.07  -5.00  12.00  29.79 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 51.93333   12.97904   4.001  0.00394 ** 
## FERTILIZER   0.81139    0.08367   9.697 1.07e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19 on 8 degrees of freedom
## Multiple R-squared:  0.9216, Adjusted R-squared:  0.9118 
## F-statistic: 94.04 on 1 and 8 DF,  p-value: 1.067e-05

summary(lm(YIELD ~ scale(FERTILIZER, scale=FALSE), data=fert))

## 
## Call:
## lm(formula = YIELD ~ scale(FERTILIZER, scale = FALSE), data = fert)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -22.79 -11.07  -5.00  12.00  29.79 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      163.50000    6.00813  27.213 3.58e-09 ***
## scale(FERTILIZER, scale = FALSE)   0.81139    0.08367   9.697 1.07e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19 on 8 degrees of freedom
## Multiple R-squared:  0.9216, Adjusted R-squared:  0.9118 
## F-statistic: 94.04 on 1 and 8 DF,  p-value: 1.067e-05

rstanarm

Using default priors

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

fert.rstanarm = stan_glm(YIELD ~ FERTILIZER, data=fert,
                         iter = 5000, warmup = 1000,
                         chains = 3, thin = 5, refresh = 0)

In the above:

• the formula and data arguments should be familiar as they are the same as for many models in R
• iter: indicates the number of MCMC iterations to perform per chain
• warmup: indicates how many of the initial MCMC samples to exclude. During the warmup stage, the MCMC sampler is tuning its sampling operations to allow it to determine the best way to create the chain. While it is in this discovery phase, the estimates it produces might be unrepresentative.
• chains: indicates the number of independent chains to run. The more chains, the more random starts. The point of multiple random starts is to mimimise the potential for the starting point to dictate the outcome. It is possible that a chain might get stuck sampling a single feature of the posterior and not discover additional features. Additional random starting points should hopefully elivate this. Either 3 or 4 is typical.
• thin: indicates the rate at which the MCMC samples should be thinned in order to eliminate autocorrelation between MCMC samples. When the step lengths between MCMC samples are small, the resulting parameter estimates will be autocorrelated (since samples close by in the chain will be more similar than those far appart). The rate of thinning necessary will depend on the degree of autocorrelation. For stan models, a thinning rate of 5 is a good starting point.
• refresh: indicates how often the terminal display should be updated. When running a large model, this can be useful as it provides a form of progress. However, when the routine is within an Rmardown document, it just results in the output including each line of progress and a lot of space is taken up unnecessarily

Having allowed rstanarm to formulate its own weakly informative priors, it is a good idea to explore what they are. Firstly, out of curiosity, it might be interesting to see what it has chosen. However, more importantly, we need to be able to document what the priors were and the rstanarm development team make it very clear that there is no guarantee that the default priors will remain the same into the future.

prior_summary(fert.rstanarm)

## Priors for model 'fert.rstanarm' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 164, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 164, scale = 160)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 2.1)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.016)
## ------
## See help('prior_summary.stanreg') for more details

This tells us:

• for the intercept, it is using a normal prior with a mean of 164 and a standard deviation of 2.5. This mean value is based on the mean of the response variable (rounded up). The 2.5 is used for all intercepts. The scale of the prior is then adjusted proportional to the variability in the response, by multiplying 2.5 by the standard deviation of the response (and rounding up).

mean(fert$YIELD)

## [1] 163.5

sd(fert$YIELD)

## [1] 63.97439

2.5*sd(fert$YIELD)

## [1] 159.936

• for the coefficients (in this case, just the slope), the default prior is a normal prior centered around 0 with a standard deviation of 2.5. This is then adjusted for the scale of the data by multiplying the 2.5 by the ratio of the standard deviation of the response to the standard deviation of the predictor (then rounded).

2.5*sd(fert$YIELD)/sd(fert$FERTILIZER)

## [1] 2.113004

• the auxillary prior represents whatever additional prior is required for the nominated model. In the case of a Gaussian model, this will be \(\sigma\), for negative binomial, it will be the reciprocal of dispersion, for gamma, it will be shape, etc . By default in rstanarm, this is a exponential with a rate of 1 which is then adjusted by devision with the standard deviation of the response.

1/sd(fert$FERTILIZER)

## [1] 0.01321157

Assessing priors

One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging preditions would ensure that the priors are unlikely to influence the actual preditions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of paramter space that are not supported by the actual data, the sampler can become unstable and have difficulty.

We can draw from the prior predictive distribution instead of conditioning on the response, by updating the model and indicating prior_PD=TRUE. After refittin the model in this way, we can plot the predictions to gain insights into the range of predictions supported by the priors alone.

fert.rstanarm1 <- update(fert.rstanarm,  prior_PD=TRUE)

ggpredict(fert.rstanarm1) %>% plot(add.data=TRUE)

## $FERTILIZER

Conclusions:

• we see that the range of predictions is faily wide and the slope could range from strongly negative to strongly positive.
• at the same time, massive yeilds are not supported, so there is some constraints on the parameters space
• nevertheless, negative yields are considered possible, even though these are clearly not logical.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

• \(\beta_0\): normal centered at 164 with a standard deviation of 10
• \(\beta_1\): normal centered at 0 with a standard deviation of 0.75
• \(\sigma\): half-cauchy centered at 0 with a standard deviation of 2

I will also overlay the raw data for comparison.

fert.rstanarm2= stan_glm(YIELD ~ FERTILIZER, data=fert,
                         prior_intercept = normal(164, 10, autoscale=FALSE),
                         prior = normal(0, 1, autoscale=FALSE),
                         prior_aux = cauchy(0, 2, autoscale=FALSE),
                         prior_PD=TRUE, 
                         iter = 5000, warmup = 1000,
                         chains = 3, thin = 5, refresh = 0
                         )

ggpredict(fert.rstanarm2) %>%
  plot(add.data=TRUE)

## $FERTILIZER

Now lets refit, conditioning on the data.

fert.rstanarm3= update(fert.rstanarm2,  prior_PD=FALSE) 

Plotting prior and posterior

posterior_vs_prior(fert.rstanarm3, color_by='vs', group_by=TRUE,
                   facet_args=list(scales='free_y'))

Conclusions:

• in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

ggpredict(fert.rstanarm3) %>% plot(add.data=TRUE)

## $FERTILIZER

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

fert.brm = brm(bf(YIELD ~ FERTILIZER), data=fert,
                         iter = 5000, warmup = 1000,
                         chains = 3, thin = 5, refresh = 0)

In the above:

• the formula and data arguments should be familiar as they are the same as for many models in R. That said, the formual specification of brms can define models that are not possible by most other routines. To facilitate this enhanced functionality, we usually define a brms formula within its own bf() function along with the family (in this case, it is gaussian, which is the default and therefore can be omitted.)
• iter: indicates the number of MCMC iterations to perform per chain
• warmup: indicates how many of the initial MCMC samples to exclude. During the warmup stage, the MCMC sampler is tuning its sampling operations to allow it to determine the best way to create the chain. While it is in this discovery phase, the estimates it produces might be unrepresentative.
• chains: indicates the number of independent chains to run. The more chains, the more random starts. The point of multiple random starts is to mimimise the potential for the starting point to dictate the outcome. It is possible that a chain might get stuck sampling a single feature of the posterior and not discover additional features. Additional random starting points should hopefully elivate this. Either 3 or 4 is typical.
• thin: indicates the rate at which the MCMC samples should be thinned in order to eliminate autocorrelation between MCMC samples. When the step lengths between MCMC samples are small, the resulting parameter estimates will be autocorrelated (since samples close by in the chain will be more similar than those far appart). The rate of thinning necessary will depend on the degree of autocorrelation. For stan models, a thinning rate of 5 is a good starting point.
• refresh: indicates how often the terminal display should be updated. When running a large model, this can be useful as it provides a form of progress. However, when the routine is within an Rmardown document, it just results in the output including each line of progress and a lot of space is taken up unnecessarily

Having allowed brms to formulate its own weakly informative priors, it is a good idea to explore what they are. Firstly, out of curiosity, it might be interesting to see what it has chosen. However, more importantly, we need to be able to document what the priors were and the brms development team make it very clear that there is no guarantee that the default priors will remain the same into the future.

prior_summary(fert.brm)

##                      prior     class       coef group resp dpar nlpar bound       source
##                     (flat)         b                                             default
##                     (flat)         b FERTILIZER                             (vectorized)
##  student_t(3, 161.5, 90.4) Intercept                                             default
##      student_t(3, 0, 90.4)     sigma                                             default

This tells us:

• for the intercept, it is using a student t (flatter normal) prior with a mean of 161.5 and a standard deviation of 90.4. These mean and standard deviation values are based on the median and median absolute deviation of the response variable (rounded).

median(fert$YIELD)

## [1] 161.5

mad(fert$YIELD)

## [1] 90.4386

standist::visualize("student_t(3,161.5,90.4)", xlim=c(-10,1000))

• for the beta coefficients (in this case, just the slope), the default prior is a inproper flat prior. A flat prior essentially means that any value between negative infinity and positive infinity are equally likely. Whilst this might seem reckless, in practice, it seems to work reasonably well for non-intercept beta parameters.

• the prior on sigma is also a student t (although only the positive half is used in the stan code), with mean of 0 and standard deviation of 90.4.

standist::visualize("student_t(3,0,90.4)", xlim=c(-10,500))

Assessing priors

One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging predictions would ensure that the priors are unlikely to influence the actual predictions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of parameter space that are not supported by the actual data, the sampler can become unstable and have difficulty.

In brms, we can inform the sampler to draw from the prior predictive distribution instead of conditioning on the response, by running the model with the sample_prior='only' argument. Unfortunately, this cannot be applied when there are flat priors (since the posteriors will necessarily extend to negative and positive infinity). Therefore, in order to use this useful routine, we need to make sure that we have defined a proper prior for all parameters.

Lets try a Gaussian (normal) distribution for the slope. We will start with what is likely to be a very wide distribution (wide with respect to what we know about the likely slope).

standist::visualize("normal(0,10)", xlim=c(-50,50))

fert.brm1 = brm(bf(YIELD ~ FERTILIZER), data=fert,
                prior=prior(normal(0, 10), class='b'), 
                sample_prior = 'only', 
                iter = 5000, warmup = 1000,
                chains = 3, thin = 5, refresh = 0)

ggpredict(fert.brm1) %>% plot(add.data=TRUE)

## $FERTILIZER

conditional_effects(fert.brm1) %>%  plot(points=TRUE)

Conclusions:

• we see that the range of predictions is faily wide and the slope could range from strongly negative to strongly positive.
• at the same time, massive yeilds are not supported, so there is some constraints on the parameters space
• nevertheless, negative yields are considered possible, even though these are clearly not logical.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

• \(\beta_0\): normal centered at 164 with a standard deviation of 10
• \(\beta_1\): normal centered at 0 with a standard deviation of 0.75
• \(\sigma\): half-cauchy centered at 0 with a standard deviation of 2
• \(\sigma\): gamma with shape parameters of 2 and 1

standist::visualize("normal(164,10)", xlim=c(-100,300))

standist::visualize("cauchy(0,2)", xlim=c(-10,10))

standist::visualize("gamma(2,1)", xlim=c(0,10))

standist::visualize("gamma(0.2,0.01)", xlim=c(0,10))

I will also overlay the raw data for comparison.

priors <- prior(normal(164, 10),  class='Intercept') +
    prior(normal(0, 1), class='b') +
    prior(gamma(2, 1),  class='sigma')

fert.brm2 = brm(bf(YIELD ~ FERTILIZER), data=fert,
                prior=priors,
                sample_prior = 'only', 
                iter = 5000, warmup = 1000,
                chains = 3, thin = 5, refresh = 0)

ggpredict(fert.brm2) %>%
  plot(add.data=TRUE)

## $FERTILIZER

conditional_effects(fert.brm2)

fert.brm3 <- update(fert.brm2,  sample_prior=TRUE, refresh=0)

Plotting prior and posterior

fert.brm3 %>% get_variables()

##  [1] "b_Intercept"     "b_FERTILIZER"    "sigma"           "prior_Intercept"
##  [5] "prior_b"         "prior_sigma"     "lp__"            "accept_stat__"  
##  [9] "stepsize__"      "treedepth__"     "n_leapfrog__"    "divergent__"    
## [13] "energy__"

fert.brm3 %>%
  posterior_samples %>%
  select(-`lp__`) %>%
  pivot_longer(cols=everything(), names_to='key', values_to='value') %>% 
  mutate(Type=ifelse(str_detect(key, 'prior'), 'Prior', 'b'),
         Class=ifelse(str_detect(key, 'Intercept'),  'Intercept',
               ifelse(str_detect(key, 'sigma'),  'Sigma',  'b'))) %>%
  ggplot(aes(x=Type,  y=value)) +
  stat_pointinterval()+
  facet_wrap(~Class,  scales='free')

Exploring the stan code

standata(fert.brm3)

## $N
## [1] 10
## 
## $Y
##  [1]  84  80  90 154 148 169 206 244 212 248
## 
## $K
## [1] 2
## 
## $X
##    Intercept FERTILIZER
## 1          1         25
## 2          1         50
## 3          1         75
## 4          1        100
## 5          1        125
## 6          1        150
## 7          1        175
## 8          1        200
## 9          1        225
## 10         1        250
## attr(,"assign")
## [1] 0 1
## 
## $prior_only
## [1] 0
## 
## attr(,"class")
## [1] "standata"

stancode(fert.brm3)

## // generated with brms 2.14.0
## functions {
## }
## data {
##   int<lower=1> N;  // total number of observations
##   vector[N] Y;  // response variable
##   int<lower=1> K;  // number of population-level effects
##   matrix[N, K] X;  // population-level design matrix
##   int prior_only;  // should the likelihood be ignored?
## }
## transformed data {
##   int Kc = K - 1;
##   matrix[N, Kc] Xc;  // centered version of X without an intercept
##   vector[Kc] means_X;  // column means of X before centering
##   for (i in 2:K) {
##     means_X[i - 1] = mean(X[, i]);
##     Xc[, i - 1] = X[, i] - means_X[i - 1];
##   }
## }
## parameters {
##   vector[Kc] b;  // population-level effects
##   real Intercept;  // temporary intercept for centered predictors
##   real<lower=0> sigma;  // residual SD
## }
## transformed parameters {
## }
## model {
##   // likelihood including all constants
##   if (!prior_only) {
##     target += normal_id_glm_lpdf(Y | Xc, Intercept, b, sigma);
##   }
##   // priors including all constants
##   target += normal_lpdf(b | 0, 0.75);
##   target += normal_lpdf(Intercept | 164, 10);
##   target += cauchy_lpdf(sigma | 0, 2)
##     - 1 * cauchy_lccdf(0 | 0, 2);
## }
## generated quantities {
##   // actual population-level intercept
##   real b_Intercept = Intercept - dot_product(means_X, b);
##   // additionally draw samples from priors
##   real prior_b = normal_rng(0,0.75);
##   real prior_Intercept = normal_rng(164,10);
##   real prior_sigma = cauchy_rng(0,2);
##   // use rejection sampling for truncated priors
##   while (prior_sigma < 0) {
##     prior_sigma = cauchy_rng(0,2);
##   }
## }

MCMC sampling diagnostics

MCMC sampling behaviour

available_mcmc()

Package	Description	function	stanarm	BRMS
Bayesplot	Traceplot	mcmc_trace	plot(mod, plotfun='trace')	mcmc_plot(mod, type='trace')
	Density plot	mcmc_dens	plot(mod, plotfun='dens')	mcmc_plot(mod, type='dens')
	Density & Trace	mcmc_combo	plot(mod, plotfun='combo')	mcmc_plot(mod, type='combo')
	ACF	mcmc_acf_bar	plot(mod, plotfun='acf_bar')	mcmc_plot(mod, type='acf_bar')
	Rhat hist	mcmc_rhat_hist	plot(mod, plotfun='rhat_hist')	mcmc_plot(mod, type='rhat_hist')
	No. Effective	mcmc_neff_hist	plot(mod, plotfun='neff_hist')	mcmc_plot(mod, type='neff_hist')
rstan	Traceplot	stan_trace	stan_trace(mod)	stan_trace(mod)
	ACF	stan_ac	stan_ac(mod)	stan_ac(mod)
	Rhat	stan_rhat	stan_rhat(mod)	stan_rhat(mod)
	No. Effective	stan_ess	stan_ess(mod)	stan_ess(mod)
	Density plot	stan_dens	stan_dens(mod)	stan_dens(mod)
ggmcmc	Traceplot	ggs_traceplot	ggs_traceplot(ggs(mod))	ggs_traceplot(ggs(mod))
	ACF	ggs_autocorrelation	ggs_autocorrelation(ggs(mod))	ggs_autocorrelation(ggs(mod))
	Rhat	ggs_Rhat	ggs_Rhat(ggs(mod))	ggs_Rhat(ggs(mod))
	No. Effective	ggs_effective	ggs_effective(ggs(mod))	ggs_effective(ggs(mod))
	Cross correlation	ggs_crosscorrelation	ggs_crosscorrelation(ggs(mod))	ggs_crosscorrelation(ggs(mod))
	Scale reduction	ggs_grb	ggs_grb(ggs(mod))	ggs_grb(ggs(mod))

In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.

There are a wide variety of tests that range from the big picture, overal chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.

rstanarm

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_data
##   mcmc_areas_ridges
##   mcmc_areas_ridges_data
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_chains_data
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_intervals_data
##   mcmc_neff
##   mcmc_neff_data
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_parcoord_data
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_data
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_data
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

• mcmc_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

plot(fert.rstanarm3, plotfun='mcmc_trace')

The chains appear well mixed and very similar

• acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

plot(fert.rstanarm3, 'acf_bar')

There is no evidence of autocorrelation in the MCMC samples

• Rhat: Rhat is a measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

plot(fert.rstanarm3, 'rhat_hist')

All Rhat values are below 1.05, suggesting the chains have converged.

• neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

plot(fert.rstanarm3, 'neff_hist')

Ratios all very high.

More diagnostics

plot(fert.rstanarm3, 'combo')

plot(fert.rstanarm3, 'violin')

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

stan_trace(fert.rstanarm3)

The chains appear well mixed and very similar

• stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

stan_ac(fert.rstanarm3) 

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

stan_rhat(fert.rstanarm3) 

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

stan_ess(fert.rstanarm3)

Ratios all very high.

stan_dens(fert.rstanarm3, separate_chains = TRUE)

ggmcmc

The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

fert.ggs <- ggs(fert.rstanarm3)
ggs_traceplot(fert.ggs)

The chains appear well mixed and very similar

• gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(fert.ggs)

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(fert.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

ggs_effective(fert.ggs)

Ratios all very high.

More diagnostics

ggs_crosscorrelation(fert.ggs)

ggs_grb(fert.ggs)

brms

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_data
##   mcmc_areas_ridges
##   mcmc_areas_ridges_data
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_chains_data
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_intervals_data
##   mcmc_neff
##   mcmc_neff_data
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_parcoord_data
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_data
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_data
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

• trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

mcmc_plot(fert.brm3, type='combo')

mcmc_plot(fert.brm3, type='trace')

mcmc_plot(fert.brm3, type='dens_overlay')

The chains appear well mixed and very similar

• acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

mcmc_plot(fert.brm3, type='acf_bar')

There is no evidence of autocorrelation in the MCMC samples

• rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

mcmc_plot(fert.brm3, type='rhat_hist')

All Rhat values are below 1.05, suggesting the chains have converged.

• neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

mcmc_plot(fert.brm3, type='neff_hist')

Ratios all very high.

More diagnostics

mcmc_plot(fert.brm3, type='combo')

mcmc_plot(fert.brm3, type='violin')

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

stan_trace(fert.brm3$fit)

The chains appear well mixed and very similar

• stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

stan_ac(fert.brm3$fit) 

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

stan_rhat(fert.brm3$fit) 

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

stan_ess(fert.brm3$fit)

Ratios all very high.

stan_dens(fert.brm3$fit, separate_chains = TRUE)

ggmcmc

The ggmcmc package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

fert.ggs <- ggs(fert.brm3, burnin=FALSE)  # does not seem to ignore burnin?
ggs_traceplot(fert.ggs)

The chains appear well mixed and very similar

• gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(fert.ggs)

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(fert.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

ggs_effective(fert.ggs)

Ratios all very high.

More diagnostics

ggs_crosscorrelation(fert.ggs)

ggs_grb(fert.ggs)

Model validation

Posterior probabilty checks

available_ppc()

Package	Description	function	stanarm	BRMS
Bayesplot	Density overlay	ppc_dens_overlay	pp_check(mod, plotfun='dens_overlay')	pp_check(mod, type='dens_overlay')
	Obs vs Pred error	ppc_error_scatter_avg	pp_check(mod, plotfun='error_scatter_avg')	pp_check(mod, type='error_scatter_avg')
	Pred error vs x	ppc_error_scatter_avg_vs_x	pp_check(mod, x=, plotfun='error_scatter_avg_vs_x')	pp_check(mod, x=, type='error_scatter_avg_vs_x')
	Preds vs x	ppc_intervals	pp_check(mod, x=, plotfun='intervals')	pp_check(mod, x=, type='intervals')
	Parial plot	ppc_ribbon	pp_check(mod, x=, plotfun='ribbon')	pp_check(mod, x=, type='ribbon')

rstanarm

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_dens_overlay_grouped
##   ppc_ecdf_overlay
##   ppc_ecdf_overlay_grouped
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_km_overlay
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_data
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

• dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

pp_check(fert.rstanarm3,  plotfun='dens_overlay')

The model draws appear to be consistent with the observed data.

• error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot

pp_check(fert.rstanarm3, plotfun='error_scatter_avg')

There is no obvious pattern in the residuals.

• error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such.

pp_check(fert.rstanarm3, x=fert$FERTILIZER, plotfun='error_scatter_avg_vs_x')

• intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

pp_check(fert.rstanarm3, x=fert$FERTILIZER, plotfun='intervals')

• ribbon: this is just an alternative way of expressing the above plot.

pp_check(fert.rstanarm3, x=fert$FERTILIZER,plotfun='ribbon')

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

#library(shinystan)
#launch_shinystan(fert.rstanarm3)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

• simulated (predicted) responses associated with each observation.
• observed values
• fitted (predicted) responses (averaged) associated with each observation

preds <- posterior_predict(fert.rstanarm3,  nsamples=250,  summary=FALSE)
fert.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = fert$YIELD,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = 'gaussian')
plot(fert.resids)

Conclusions:

• the simulated residuals do not suggest any issues with the residuals
• there is no evidence of a lack of fit.

brms

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_dens_overlay_grouped
##   ppc_ecdf_overlay
##   ppc_ecdf_overlay_grouped
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_km_overlay
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_data
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

• dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

pp_check(fert.brm3,  type='dens_overlay')

The model draws appear to be consistent with the observed data.

• error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot

pp_check(fert.brm3, type='error_scatter_avg')

There is no obvious pattern in the residuals.

• error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such.

pp_check(fert.brm3, x='FERTILIZER', type='error_scatter_avg_vs_x')

• intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

pp_check(fert.brm3, x='FERTILIZER', type='intervals')

• ribbon: this is just an alternative way of expressing the above plot.

pp_check(fert.brm3, x='FERTILIZER',type='ribbon')

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

#library(shinystan)
#launch_shinystan(fert.brm3)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

• simulated (predicted) responses associated with each observation.
• observed values
• fitted (predicted) responses (averaged) associated with each observation

preds <- posterior_predict(fert.brm3,  nsamples=250,  summary=FALSE)
fert.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = fert$YIELD,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = 'gaussian')
plot(fert.resids)

Conclusions:

• the simulated residuals do not suggest any issues with the residuals
• there is no evidence of a lack of fit.

Partial effects plots

rstanarm

ggpredict

fert.rstanarm3 %>% ggpredict() %>% plot(add.data=TRUE)

## $FERTILIZER

ggemmeans

fert.rstanarm3 %>% ggemmeans(~FERTILIZER) %>% plot(add.data=TRUE)

fitted_draws

fert.rstanarm3 %>% fitted_draws(newdata=fert) %>%
 median_hdci() %>%
 ggplot(aes(x=FERTILIZER, y=.value)) +
 geom_ribbon(aes(ymin=.lower, ymax=.upper), fill='blue', alpha=0.3) + 
 geom_line()

brms

conditional_effects

fert.brm3 %>% conditional_effects() 

fert.brm3 %>% conditional_effects(spaghetti=TRUE,nsamples=200) 

ggpredict

fert.brm3 %>% ggpredict() %>% plot(add.data=TRUE)

## $FERTILIZER

ggemmeans

fert.brm3 %>% ggemmeans(~FERTILIZER) %>% plot(add.data=TRUE)

fitted_draws

fert.brm3 %>%
    fitted_draws(newdata=fert) %>%
    median_hdci() %>%
    ggplot(aes(x=FERTILIZER, y=.value)) +
    geom_ribbon(aes(ymin=.lower, ymax=.upper), fill='blue', alpha=0.3) + 
    geom_line()

Model investigation

Package	Function	Description
	as.matrix()	Returns \(n\times p\) matrix
tidybayes	tidy_draws()	Returns \(n\times p\) tibble with addition info about the chain, iteration and draw
tidybayes	spread_draws()	Returns \(n\times r\) tibble (where \(r\) is the number of requested parameters) with additional info about chain, iteraction and draw
tidybayes	gather_draws()	Returns a gathered spread_draws tibble with additional info about chain, iteraction and draw
brms	posterior_samples()	Returns \(n\times p\) data.frame

where \(n\) is the number of MCMC samples and \(p\) is the number of parameters to estimate.

Function	Description
median_qi	Median and quantiles
median_hdi	Median and Highest Probability Density Interval
median_hdci	Median and continuous Highest Probability Density Interval

It is important to remember that when working with Bayesian models, everything is a distribution. Whilst point estimates (such as a mean) of the parameters can be calculated from these distributions, we start off with a large number of estimates per parameter.

rstanarm

rstanarm captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

summary(fert.rstanarm3)

## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      YIELD ~ FERTILIZER
##  algorithm:    sampling
##  sample:       2400 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 10
##  predictors:   2
## 
## Estimates:
##               mean   sd   10%   50%   90%
## (Intercept) 52.6   13.0 36.6  52.8  68.3 
## FERTILIZER   0.8    0.1  0.7   0.8   0.9 
## sigma       19.4    5.1 13.9  18.4  26.2 
## 
## Fit Diagnostics:
##            mean   sd    10%   50%   90%
## mean_PPD 163.6    8.4 153.3 163.5 173.8
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.3  1.0  2242 
## FERTILIZER    0.0  1.0  2442 
## sigma         0.1  1.0  2278 
## mean_PPD      0.2  1.0  2426 
## log-posterior 0.0  1.0  2181 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Conclusions:

• in the Model Info, we are informed that the total MCMC posterior sample size is 2400 and that there were 10 raw observations.
• the estimated intercept (expected yield when fertilizer concentration is 0) is 52.59. This is the mean of the posterior distribution for this parameter.
• the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.81 (mean) or 0.7 (median) with a standard deviation of 0. The 90% credibility intervals indicate that we are 90% confident that the slope is between 0.81 and 0.8 - e.g. there is a significant positive trend.
• sigma is estimated to be 19.42
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

tidyMCMC

tidyMCMC(fert.rstanarm3$stanfit, estimate.method='median',  conf.int=TRUE,
         conf.method='HPDinterval',  rhat=TRUE, ess=TRUE)

Conclusions:

• the estimated intercept (expected yield when fertilizer concentration is 0) is 52.82. This is the median of the posterior distribution for this parameter.
• the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.8 (median) with a standard error of 0.09. The 95% credibility intervals indicate that we are 95% confident that the slope is between 0.63 and 0.98 - e.g. there is a significant positive trend.
• sigma is estimated to be 18.43
• Rhat and number of effective samples for each parameter are also provided and all look good.

gather_draws

fert.draw <- fert.rstanarm3 %>% gather_draws(`(Intercept)`, FERTILIZER, sigma)
## OR via regex
fert.draw <- fert.rstanarm3 %>% gather_draws(`.Intercept.*|FERT.*|sigma`,  regex=TRUE)
fert.draw

We can then summarise this

fert.draw %>% median_hdci

fert.rstanarm3 %>% 
  gather_draws(`(Intercept)`, FERTILIZER, sigma) %>% 
  ggplot() + 
  stat_halfeye(aes(x=.value,  y=.variable)) +
  facet_wrap(~.variable, scales='free')

Conclusions:

• the estimated intercept (expected yield when fertilizer concentration is 0) is 52.82. This is the median of the posterior distribution for this parameter.
• the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.8 (mean) or (median). The 95% credibility intervals indicate that we are 95% confident that the slope is between 0.63 and 0.98 - e.g. there is a significant positive trend.
• sigma is estimated to be 18.43

bayesplot

fert.rstanarm3 %>% plot(plotfun='mcmc_intervals') 

tidy_draws

This is purely a graphical depiction on the posteriors.

fert.rstanarm3 %>% tidy_draws()

spread_draws

fert.rstanarm3 %>% spread_draws(`(Intercept)`, FERTILIZER, sigma)

# OR via regex
fert.rstanarm3 %>% spread_draws(`.Intercept.*|FERT.*|sigma`,  regex=TRUE)

posterior_samples

fert.rstanarm3 %>% posterior_samples() %>% as_tibble()

\(R^2\)

fert.rstanarm3 %>% bayes_R2() %>% median_hdci

Empirical p-value

It is possible to obtain an empirical p-value (for those that cannot live without p-values). This is essentially compares (for any posterior) that the mean is zero compared to a multivatiate distribution with elliptical contours.

mcmcpvalue <- function(samp)
{
    ## elementary version that creates an empirical p-value for the
    ## hypothesis that the columns of samp have mean zero versus a
    ## general multivariate distribution with elliptical contours.
    
    ## differences from the mean standardized by the observed
    ## variance-covariance factor
    
    ## Note, I put in the bit for single terms
    if (length(dim(samp))==0) {
        std <- backsolve(chol(var(samp)),cbind(0, t(samp)) - mean(samp),transpose = TRUE)
        sqdist <- colSums(std * std)
        sum(sqdist[-1] > sqdist[1])/length(samp)
    }
    else {
        std <- backsolve(chol(var(samp)),cbind(0, t(samp)) - colMeans(samp),transpose = TRUE)
        sqdist <- colSums(std * std)
        sum(sqdist[-1] > sqdist[1])/nrow(samp)
    }
    
}

mcmcpvalue(as.matrix(fert.rstanarm3)[, c("FERTILIZER")])

## [1] 0

brms

brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

summary(fert.brm3)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: YIELD ~ FERTILIZER 
##    Data: fert (Number of observations: 10) 
## Samples: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup samples = 2400
## 
## Population-Level Effects: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     53.35     12.89    28.13    78.88 1.00     2395     2274
## FERTILIZER     0.80      0.09     0.62     0.96 1.00     2415     2270
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    19.21      5.10    12.18    31.96 1.00     2266     2239
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

• in the initial block of information, we are reminded on the formula as well as chain fitting characteristics. We are also informed that the total number of post-warmup MCMC samples is is 2400.
• the estimated intercept (expected yield when fertilizer concentration is 0) is 53.35. This is the mean of the posterior distribution for this parameter.
• the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.8 (mean) with a standard error of 0.09. The 95% credibility intervals indicate that we are 95% confident that the slope is between 0.62 and 0.96 - e.g. there is a significant positive trend.
• sigma is estimated to be 19.21
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

tidyMCMC

tidyMCMC(fert.brm3$fit, estimate.method='median',  conf.int=TRUE,  conf.method='HPDinterval',  rhat=TRUE, ess=TRUE)

Conclusions:

• the estimated intercept (expected yield when fertilizer concentration is 0) is 53.15. This is the median of the posterior distribution for this parameter.
• the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 0.8 (median) with a standard error of 0.09. The 95% credibility intervals indicate that we are 95% confident that the slope is between 0.64 and 0.97 - e.g. there is a significant positive trend.
• sigma is estimated to be 18.25
• Rhat and number of effective samples for each parameter are also provided and all look good.

gather_draws

fert.brm3 %>% get_variables()

##  [1] "b_Intercept"     "b_FERTILIZER"    "sigma"           "prior_Intercept"
##  [5] "prior_b"         "prior_sigma"     "lp__"            "accept_stat__"  
##  [9] "stepsize__"      "treedepth__"     "n_leapfrog__"    "divergent__"    
## [13] "energy__"

fert.brm3 %>%
    gather_draws(b_Intercept, b_FERTILIZER, sigma) %>%
    median_hdci

## OR via regex
fert.brm3 %>%
    gather_draws(`b_.*|sigma`,  regex=TRUE) %>%
    median_hdci

fert.brm3 %>% 
  gather_draws(b_Intercept, b_FERTILIZER, sigma) %>% 
  ggplot() + 
  stat_halfeye(aes(x=.value,  y=.variable)) +
  facet_wrap(~.variable, scales='free')

fert.brm3 %>% 
  gather_draws(b_Intercept, b_FERTILIZER, sigma) %>% 
  ggplot() + 
  stat_halfeye(aes(x=.value,  y=.variable)) 

Conclusions:

• the estimated intercept (expected yield when fertilizer concentration is 0) is 0.8. This is the median of the posterior distribution for this parameter.
• the estimated slope (rate at which yield changes per 1 unit change in fertilizer concentration), is 53.15 (mean) or (median). The 95% credibility intervals indicate that we are 95% confident that the slope is between 27.1 and 77.24 - e.g. there is a significant positive trend.
• sigma is estimated to be 18.25

bayesplot

fert.brm3 %>% mcmc_plot(type='intervals') 

tidy_draws

Return the draws (samples) for each parameter

fert.brm3 %>% tidy_draws()

spread_draws

Select just the parameters of interest.

fert.brm3 %>% spread_draws(b_Intercept, b_FERTILIZER, sigma)

# OR via regex
fert.brm3 %>% spread_draws(`b_.*|sigma`,  regex=TRUE)

posterior_samples

Just the main parameters

fert.brm3 %>% posterior_samples() %>% as_tibble()

\(R^2\)

fert.brm3 %>% bayes_R2(summary=FALSE) %>% median_hdci

Empirical p-value

It is possible to obtain an empirical p-value (for those that cannot live without p-values). This is essentially compares (for any posterior) that the mean is zero compared to a multivatiate distribution with elliptical contours.

mcmcpvalue <- function(samp)
{
    ## elementary version that creates an empirical p-value for the
    ## hypothesis that the columns of samp have mean zero versus a
    ## general multivariate distribution with elliptical contours.
    
    ## differences from the mean standardized by the observed
    ## variance-covariance factor
    
    ## Note, I put in the bit for single terms
    if (length(dim(samp))==0) {
        std <- backsolve(chol(var(samp)),cbind(0, t(samp)) - mean(samp),transpose = TRUE)
        sqdist <- colSums(std * std)
        sum(sqdist[-1] > sqdist[1])/length(samp)
    }
    else {
        std <- backsolve(chol(var(samp)),cbind(0, t(samp)) - colMeans(samp),transpose = TRUE)
        sqdist <- colSums(std * std)
        sum(sqdist[-1] > sqdist[1])/nrow(samp)
    }
    
}

mcmcpvalue(as.matrix(fert.brm3)[, c("b_FERTILIZER")])

## [1] 0

Predictions

There are a large number of condidate routines for performing prediction. We will go through many of these. It is worth noting that prediction is technically the act of estimating what we expect to get if we were to collect a single new observation from a particular population (e.g. a specific level of fertilizer concentration). Often this is not what we want. Often we want the fitted values - estimates of what we expect to get if we were to collect multiple new observations and average them.

So while fitted values represent the expected underlying processes occuring in the system, predicted values represent our expectations from sampling from such processes.

Package	Function	Description	Summarise with
rstantools	posterior_predict	Draw from the posterior of a prediction (includes sigma) - predicts single observations	tidyMCMC()
rstantools	posterior_linpred	Draw from the posterior of the fitted values (on the link scale) - predicts average observations	tidyMCMC()
rstantools	posterior_epred	Draw from the posterior of the fitted values (on the response scale) - predicts average observations	tidyMCMC()
tidybayes	add_predicted_draws	Adds draws from the posterior of predictions to a data frame (of prediction data)	median_hdci()
tidybayes	add_fitted_draws	Adds draws from the posterior of fitted values to a data frame (of prediction data)	median_hdci()

For simple models prediction is essentially taking the model formula complete with parameter (coefficient) estimates and solving for new values of the predictor. To explore this, we will use the fitted model to predict Yield for a Fertilizer concentration of 110.

We will therefore start by establishing this prediction domain as a data frame to use across all of the prediction routines.

## establish a data set that defines the new data to predict against
newdata = data.frame(FERTILIZER=110)

rstanarm

fitted_draws

fert.rstanarm3 %>% fitted_draws(newdata=newdata) %>%
  median_hdci

emmeans

fert.rstanarm3 %>% emmeans(~FERTILIZER,  at=newdata)

##  FERTILIZER emmean lower.HPD upper.HPD
##         110    141       130       152
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

posterior_predict

fert.rstanarm3 %>% posterior_predict(newdata=newdata) %>%
  tidyMCMC(estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

posterior_epred

fert.rstanarm3 %>% posterior_epred(newdata=newdata) %>%
  tidyMCMC(estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

posterior_linpred

fert.rstanarm3 %>% posterior_linpred(newdata=newdata) %>%
  tidyMCMC(estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

add_predicted_draws

newdata %>% add_predicted_draws(fert.rstanarm3) %>%
  median_hdci

add_fitted_draws

newdata %>% add_fitted_draws(fert.rstanarm3) %>%
  median_hdci

manual calculations

## Establish an appropriate model matrix
Xmat = model.matrix(~FERTILIZER, data=newdata)
### get the posterior draws for the linear predictor
coefs <- fert.rstanarm3 %>% posterior_samples(pars=c('(Intercept)','FERTILIZER')) %>% as.matrix()
fit <- coefs %*% t(Xmat)
fit %>% median_hdci

brms

fitted_draws

fert.brm3 %>% fitted_draws(newdata=newdata) %>%
  median_hdci

emmeans

fert.brm3 %>% emmeans(~FERTILIZER,  at=newdata)

##  FERTILIZER emmean lower.HPD upper.HPD
##         110    142       131       153
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

posterior_predict

fert.brm3 %>% posterior_predict(newdata=newdata) %>% as.mcmc() %>%
  tidyMCMC(estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

posterior_epred

fert.brm3 %>% posterior_epred(newdata=newdata) %>% as.mcmc() %>%
  tidyMCMC(estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

posterior_linpred

fert.brm3 %>% posterior_linpred(newdata=newdata) %>% as.mcmc() %>%
  tidyMCMC(estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

add_predicted_draws

newdata %>% add_predicted_draws(fert.brm3) %>%
  median_hdci

add_fitted_draws

newdata %>% add_fitted_draws(fert.brm3) %>%
  median_hdci

manual calculations

## Establish an appropriate model matrix
Xmat = model.matrix(~FERTILIZER, data=newdata)
### get the posterior draws for the linear predictor
coefs <- fert.brm3 %>% posterior_samples(pars=c('b_Intercept','b_FERTILIZER')) %>% as.matrix()
fit <- coefs %*% t(Xmat)
fit %>% median_hdci

Hypothesis testing

Since we have the entire posterior, we are able to make probablity statements. We simply count up the number of MCMC sample draws that satisfy a condition (e.g represent a slope greater than 0) and then divide by the total number of MCMC samples.

For this exersize, we will explore the following:

• what this the probability that Yield increases with increasing Fertilizer concentration - this is just asking what proportion of the drawn slopes are greater than 0.
• how much of a change in yield do we expect if fertilizer concentration is increased from 100 to 200 units and what is the probability that the change is more than 50%

rstanarm

Prob. of effect

fert.rstanarm3 %>% hypothesis('FERTILIZER>0')

## Hypothesis Tests for class :
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (FERTILIZER) > 0     0.81      0.09     0.67     0.95        Inf         1
##   Star
## 1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Conclusions:

• the probability of a positive effect of Fertilizer on grass Yield is 1. That is, there is, we are 100% confident that there is an effect.
• the evidence ratio is normally the ratio of the number of cases that satisty the hypothesis to the number of cases that do not. Since the number of cases that do not satisfy the hypothesis is 0, the evidence ratio is Inf (since division by 0)

Alternatively…

fert.rstanarm3 %>% tidy_draws() %>% summarise(P=sum(FERTILIZER>0)/n())

## # A tibble: 1 x 1
##       P
##   <dbl>
## 1     1

Conclusions:

• the probability of a positive effect of Fertilizer on grass Yield is 1. That is, there is, we are 100% confident that there is an effect.

Change 100 to 200

newdata = list(FERTILIZER=c(200, 100)) 
fert.rstanarm3 %>% emmeans(~FERTILIZER,  at=newdata) %>% pairs()

##  contrast  estimate lower.HPD upper.HPD
##  200 - 100     80.5      63.3      97.8
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

Conclusions:

• the change in Yield associated with an increase in Fertilizer concentration from 100 to 200 is 80.49.

If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.

fert.mcmc <- fert.rstanarm3 %>% emmeans(~FERTILIZER,  at=newdata) %>% 
  tidy_draws() %>%
  rename_with(~str_replace(., 'FERTILIZER ', 'p')) %>%
  mutate(Eff=p200 - p100,
         PEff=100*Eff/p100)
fert.mcmc %>% head

Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).

fert.mcmc %>% tidyMCMC(estimate.method='median',
                       conf.int=TRUE, conf.method='HPDinterval')

Alternatively, we could use median_hdci

fert.mcmc %>% median_hdci(PEff)

Conclusions:

• the percentage change in Yield associated with an increase in Fertilizer concentration from 100 to 200 is 60.4% with a 95% credibility interval of 44.37 76.5

To get the probability that the effect is greater than a 50% increase.

fert.mcmc %>% summarise(P=sum(PEff>50)/n())

## # A tibble: 1 x 1
##       P
##   <dbl>
## 1 0.918

Conclusions:

• the probability that Yield will increase by more than 50 following an increase in Fertilizer concentration from 100 to 200 is 0.92.

Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).

fert.mcmc %>% hypothesis('PEff>50')

## Hypothesis Tests for class :
##        Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (PEff)-(50) > 0    10.83      8.37    -2.24     24.5      11.18      0.92
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

brms

Prob. of effect

fert.brm3 %>% hypothesis('FERTILIZER>0')

## Hypothesis Tests for class b:
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (FERTILIZER) > 0      0.8      0.09     0.66     0.94        Inf         1
##   Star
## 1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Conclusions:

• the probability of a positive effect of Fertilizer on grass Yield is 1. That is, there is, we are 100% confident that there is an effect.
• the evidence ratio is normally the ratio of the number of cases that satisty the hypothesis to the number of cases that do not. Since the number of cases that do not satisfy the hypothesis is 0, the evidence ratio is Inf (since division by 0)

Alternatively…

fert.brm3 %>% tidy_draws() %>% summarise(P=sum(b_FERTILIZER>0)/n())

## # A tibble: 1 x 1
##       P
##   <dbl>
## 1     1

Conclusions:

• the probability of a positive effect of Fertilizer on grass Yield is 1. That is, there is, we are 100% confident that there is an effect.

Change 100 to 200

newdata = list(FERTILIZER=c(200, 100)) 
fert.brm3 %>% emmeans(~FERTILIZER,  at=newdata) %>% pairs()

##  contrast  estimate lower.HPD upper.HPD
##  200 - 100     80.4      63.8      97.4
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

Conclusions:

• the change in Yield associated with an increase in Fertilizer concentration from 100 to 200 is 80.44.

If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.

fert.mcmc <- fert.brm3 %>% emmeans(~FERTILIZER,  at=newdata) %>% 
  tidy_draws() %>%
  rename_with(~str_replace(., 'FERTILIZER ', 'p')) %>%
  mutate(Eff=p200 - p100,
         PEff=100*Eff/p100)
fert.mcmc %>% head

Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).

fert.mcmc %>% tidyMCMC(estimate.method='median',
                       conf.int=TRUE, conf.method='HPDinterval')

Alternatively, we could use median_hdci

fert.mcmc %>% median_hdci(PEff)

Conclusions:

• the percentage change in Yield associated with an increase in Fertilizer concentration from 100 to 200 is 60.19% with a 95% credibility interval of 45.51 77.73

To get the probability that the effect is greater than a 50% increase.

fert.mcmc %>% summarise(P=sum(PEff>50)/n())

## # A tibble: 1 x 1
##       P
##   <dbl>
## 1 0.899

Conclusions:

• the probability that Yield will increase by more than 50 following an increase in Fertilizer concentration from 100 to 200 is 0.9.

Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).

fert.mcmc %>% hypothesis('PEff>50')

## Hypothesis Tests for class :
##        Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (PEff)-(50) > 0    10.32      8.27    -2.71     23.9       8.88       0.9
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Summary figures

rstanarm

fert.list = with(fert, list(FERTILIZER = seq(min(FERTILIZER), max(FERTILIZER), len=100)))
newdata = emmeans(fert.rstanarm3, ~FERTILIZER, at=fert.list) %>% as.data.frame
head(newdata)

ggplot(newdata, aes(y=emmean, x=FERTILIZER)) + 
geom_point(data=fert, aes(y=YIELD)) +
geom_line() + 
geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD), fill='blue', alpha=0.3) +
scale_y_continuous('YIELD') +
scale_x_continuous('FERTILIZER') +
theme_classic()

## spaghetti plot
newdata = emmeans(fert.rstanarm3, ~FERTILIZER, at=fert.list) %>%
  gather_emmeans_draws()
newdata %>% head

ggplot(newdata,  aes(y=.value,  x=FERTILIZER)) +
  geom_line(aes(group=.draw),  alpha=0.01) +
  geom_point(data=fert,  aes(y=YIELD))

brms

fert.list = with(fert, list(FERTILIZER = seq(min(FERTILIZER), max(FERTILIZER), len=100)))
newdata = emmeans(fert.brm3, ~FERTILIZER, at=fert.list) %>% as.data.frame
head(newdata)

ggplot(newdata, aes(y=emmean, x=FERTILIZER)) + 
geom_point(data=fert, aes(y=YIELD)) +
geom_line() + 
geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD), fill='blue', alpha=0.3) +
scale_y_continuous('YIELD') +
scale_x_continuous('FERTILIZER') +
theme_classic()

## spaghetti plot
newdata = emmeans(fert.brm3, ~FERTILIZER, at=fert.list) %>%
  gather_emmeans_draws()
newdata %>% head

ggplot(newdata,  aes(y=.value,  x=FERTILIZER)) +
  geom_line(aes(group=.draw),  alpha=0.01) +
  geom_point(data=fert,  aes(y=YIELD))

References

Fowler, J., L. Cohen, and P. Jarvis. 1998. Practical Statistics for Field Biology. England: John Wiley & Sons.